#include "qlib/geo/path.h"

#include <gdstk/gdstk.hpp>

#include <cmath>

namespace qlib {
Path::Path() {
  // memset((void*)this, 0, sizeof(Path));

  simple_path = false;
  scale_width = false;
  owner = nullptr;

  spine.tolerance = 0.001;
  width = 0;
  radius = 0;

  elements.resize(1);
  auto& ele = elements[0];

  ele.half_width_and_offset.resize(1);
  ele.half_width_and_offset[0] = qlib::Vec2d(0, 0);

  elements[0].bend_type = BendType::Circular;
  elements[0].join_type = JoinType::Round;

  type = PathType::Spine;
}

Path::Path(const Path& oth) {
  operator=(oth);
}
Path& Path::operator=(const Path& other) {
  type = other.type;
  spine = other.spine;
  box = other.box;
  width = other.width;
  center_width = other.center_width;
  end_type = other.end_type;
  end_extensions = other.end_extensions;
  layer = other.layer;
  datatype = other.datatype;
  elements = other.elements;
  simple_path = other.simple_path;
  scale_width = other.scale_width;
  m_raith_data = other.m_raith_data;
  repetition = other.repetition;
  m_net_name = other.m_net_name;
  radius = other.radius;

  properties.resize(other.properties.size());
  for (size_t i = 0; i < properties.size(); i++) {
    if (other.properties[i]) {
      properties[i] = Property::create();
      *properties[i] = *other.properties[i];
    }
  }

  polys = other.polys;

  return *this;
}

Path::s_ptr Path::copy() const {
  auto res = Path::create();
  *res = *this;

  return res;
}

void Path::init(const Vec2d initial_position, double width, double offset, double tolerance, uint32_t _layer, uint32_t _datatype) {
  spine.tolerance = tolerance;
  spine.append(initial_position);
  width /= 2;

  layer = _layer;
  datatype = _datatype;
  for (uint64_t i = 0; i < elements.size(); i++) {
    elements[i].half_width_and_offset.push_back(Vec2d{width, offset});
  }
}
void Path::init(const Vec2d initial_position, const double* width, const double* offset, double tolerance, uint32_t _layer, uint32_t _datatype) {
  spine.tolerance = tolerance;
  spine.append(initial_position);

  layer = _layer;
  datatype = _datatype;
  for (uint64_t i = 0; i < elements.size(); i++) {
    elements[i].half_width_and_offset.push_back(Vec2d{0.5 * width[i], offset[i]});
  }
}
void Path::init(const Vec2d initial_position, uint64_t num_elements_, double width, double separation, double tolerance, uint32_t _layer,
                uint32_t _datatype) {
  elements.resize(num_elements_);
  spine.tolerance = tolerance;
  spine.append(initial_position);
  width /= 2;
  double i0 = 0.5 * (num_elements_ - 1);

  layer = _layer;
  datatype = _datatype;
  for (uint64_t i = 0; i < elements.size(); i++) {
    elements[i].half_width_and_offset.push_back(Vec2d{width, separation * (i - i0)});
  }
}
void Path::init(const Vec2d initial_position, uint64_t num_elements_, const double* width, const double* offset, double tolerance, uint32_t _layer,
                uint32_t _datatype) {
  elements.resize(num_elements_);
  spine.tolerance = tolerance;
  spine.append(initial_position);

  layer = _layer;
  datatype = _datatype;
  for (uint64_t i = 0; i < elements.size(); i++) {
    elements[i].half_width_and_offset.push_back(Vec2d{0.5 * width[i], offset[i]});
  }
}

void Path::init_spine(Vec2dArray points, uint32_t _layer, uint32_t _datatype) {
  type = PathType::Spine;
  elements.resize(1);
  spine.tolerance = 0;
  spine.points = points;

  layer = _layer;
  datatype = _datatype;
  for (uint64_t i = 0; i < 1; i++) {
    elements[i].half_width_and_offset.push_back(Vec2d{0, 0});
  }
  updateType();
}

void Path::init_waveguide(Vec2dArray points, double line_width, double center_width, double radius, double tolerance, uint32_t _layer,
                          uint32_t _datatype) {
  type = PathType::WaveGuide;
  elements.resize(2);
  spine.tolerance = tolerance;
  spine.points = points;
  auto offset = center_width / 2 + line_width / 2;

  layer = _layer;
  datatype = _datatype;
  for (uint64_t i = 0; i < elements.size(); i++) {
    elements[i].half_width_and_offset.push_back(Vec2d{line_width, offset});
  }
  updateType();
}

void Path::apply_repetition(PathPtrArray& result) {
  if (repetition.type == RepetitionType::None)
    return;

  Array<Vec2d> offsets = {};
  repetition.get_offsets(offsets);
  repetition = Repetition();

  // Skip first offset (0, 0)
  Vec2d* offset_p = offsets.data() + 1;
  vector_reserve_now(result, offsets.size() - 1);

  for (uint64_t offset_count = offsets.size() - 1; offset_count > 0; offset_count--) {
    Path::s_ptr path = Path::create();
    *path = *this;
    path->translate(offset_p->x, offset_p->y);
    offset_p++;
    result.push_back(path);
  }

  offsets.clear();
  return;
}
void Path::arc(Vec2d radius, double initial_angle, double final_angle, double rotation, const double* width, const double* offset) {
  spine.arc(radius.x, radius.y, initial_angle, final_angle, rotation);
  fill_offsets_and_widths(width, offset);
}

void Path::bezier(const Vec2dArray& point_array, const double* width, const double* offset, bool relative) {
  spine.bezier(point_array, relative);
  fill_offsets_and_widths(width, offset);
}

void Path::fill_offsets_and_widths(const double* width, const double* offset) {
  if (elements.size() < 1)
    return;
  const uint64_t num_pts = spine.points.size() - elements[0].half_width_and_offset.size();
  for (uint64_t ne = 0; ne < elements.size(); ne++) {
    Array<Vec2d>* half_width_and_offset = &elements[ne].half_width_and_offset;
    const Vec2d initial_widoff = (*half_width_and_offset)[half_width_and_offset->size() - 1];
    const double wid = width == NULL ? 0 : 0.5 * *width++ - initial_widoff.u;
    const double off = offset == NULL ? 0 : *offset++ - initial_widoff.v;
    const Vec2d widoff_change = Vec2d{wid, off};

    vector_reserve_now(*half_width_and_offset, num_pts);
    for (uint64_t i = 1; i <= num_pts; i++)
      half_width_and_offset->push_back(initial_widoff + widoff_change * ((double)i / num_pts));
  }
}

void Path::remove_overlapping_points() {
  const double tol_sq = spine.tolerance * spine.tolerance;
  Array<Vec2d>* array = &spine.points;
  for (uint64_t i = 1; i < array->size();) {
    double dis = ((*array)[i] - (*array)[i - 1]).length_sq();
    if (dis <= tol_sq) {
      array_remove(*array, i);
      LineElement* el = elements.data();
      for (uint64_t ne = 0; ne < elements.size(); ne++, el++)
        array_remove(el->half_width_and_offset, i);
    } else {
      i++;
    }
  }
}

void Path::path_spines(Vec2dArray& res) {
  res = spine.points;
}

ErrorCode Path::element_center(const LineElement* el, Array<Vec2d>& result) {
  const Array<Vec2d> spine_points = spine.points;
  const BendType bend_type = el->bend_type;
  const double bend_radius = el->bend_radius;
  const double path_half_widths = *(double*)el->half_width_and_offset.data();
  const double path_offsets = *((double*)el->half_width_and_offset.data() + 1);
  Vec2d spine_normal = (spine_points[1] - spine_points[0]).ortho();
  spine_normal.normalize();
  Vec2d p0 = spine_points[0] + spine_normal * path_offsets;
  Vec2d p1 = spine_points[1] + spine_normal * path_offsets;
  Vec2d t0 = p1 - p0;
  t0.normalize();
  Vec2d n0 = t0.ortho();
  result.push_back(p0);

  if (spine_points.size() > 2) {
    Curve arc_curve = {};
    arc_curve.tolerance = spine.tolerance;
    spine_normal = (spine_points[2] - spine_points[1]).ortho();
    spine_normal.normalize();
    Vec2d p2 = spine_points[1] + spine_normal * path_offsets;
    Vec2d p3 = spine_points[2] + spine_normal * path_offsets;
    Vec2d t1 = p3 - p2;
    t1.normalize();
    Vec2d n1 = t1.ortho();
    double u0, u1;
    segments_intersection(p1, t0, p2, t1, u0, u1);
    Vec2d p_next = 0.5 * (p1 + u0 * t0 + p2 + u1 * t1);
    Vec2d p = p0;
    double len_next = (p_next - p).length();

    for (uint64_t i = 1; i < spine_points.size() - 1; i++) {
      Vec2d t2, n2;
      p0 = p2;
      p1 = p3;
      p = p_next;

      if (i + 2 == spine_points.size()) {
        // Last point: no need to find an intersection
        p_next = p1;
        t2 = Vec2d{0, 0};
        n2 = Vec2d{0, 0};
      } else {
        spine_normal = (spine_points[i + 2] - spine_points[i + 1]).ortho();
        spine_normal.normalize();
        p2 = spine_points[i + 1] + spine_normal * path_offsets;
        p3 = spine_points[i + 2] + spine_normal * path_offsets;
        t2 = p3 - p2;
        t2.normalize();
        n2 = t2.ortho();
        segments_intersection(p1, t1, p2, t2, u0, u1);
        p_next = 0.5 * (p1 + u0 * t1 + p2 + u1 * t2);
      }

      if (bend_type != BendType::None) {
        const double bend_dir = t0.cross(t1) < 0 ? -1 : 1;
        const double len_prev = len_next;
        len_next = (p_next - p).length();
        const Vec2d sum_t = t0 + t1;
        const double len_factor = (fabs(sum_t.x) > fabs(sum_t.y)) ? bend_dir * (n0.x - n1.x) / sum_t.x : bend_dir * (n0.y - n1.y) / sum_t.y;
        const double radius = bend_radius - bend_dir * path_offsets;
        const double len_required = len_factor * radius;
        if (len_required > len_prev || len_required > len_next || radius <= path_half_widths) {
          // Not enough room for the bend
          result.push_back(p);
        } else {
          len_next -= len_required;
          if (bend_dir < 0) {  // Right turn
            const double initial_angle = n0.angle();
            double final_angle = n1.angle();
            if (final_angle > initial_angle)
              final_angle -= 2 * M_PI;
            const Vec2d center = p - 0.5 * radius * (n0 + n1 + len_factor * (t0 - t1));
            if (bend_type == BendType::Circular) {
              const Vec2d arc_start = center + n0 * radius;
              arc_curve.append(arc_start);
              arc_curve.arc(radius, radius, initial_angle, final_angle, 0);
              array_extend(result, arc_curve.points);
              arc_curve.points.clear();
            } else if (bend_type == BendType::Function) {
              Array<Vec2d> bend_array = (*el->bend_function)(radius, initial_angle, final_angle, center, el->bend_function_data);
              array_extend(result, bend_array);
              bend_array.clear();
            }
          } else {  // Left turn
            const double initial_angle = (-n0).angle();
            double final_angle = (-n1).angle();
            if (final_angle < initial_angle)
              final_angle += 2 * M_PI;
            Vec2d center = p + 0.5 * radius * (n0 + n1 + len_factor * (t1 - t0));
            if (bend_type == BendType::Circular) {
              const Vec2d arc_start = center - n0 * radius;
              arc_curve.append(arc_start);
              arc_curve.arc(radius, radius, initial_angle, final_angle, 0);
              array_extend(result, arc_curve.points);
              arc_curve.points.clear();
            } else if (bend_type == BendType::Function) {
              Array<Vec2d> bend_array = (*el->bend_function)(radius, initial_angle, final_angle, center, el->bend_function_data);
              array_extend(result, bend_array);
              bend_array.clear();
            }
          }
        }
      } else {
        result.push_back(p);
      }

      t0 = t1;
      n0 = n1;
      t1 = t2;
      n1 = n2;
    }
    //arc_curve.clear();
  }
  result.push_back(p1);
  return ErrorCode::NoError;
}

uint64_t Path::commands(const CurveInstruction* items, uint64_t count) {
  uint64_t result = spine.commands(items, count);
  fill_offsets_and_widths(NULL, NULL);
  return result;
}

void Path::cubic(const Vec2dArray& point_array, const double* width, const double* offset, bool relative) {
  spine.cubic(point_array, relative);
  fill_offsets_and_widths(width, offset);
}
void Path::cubic_smooth(const Vec2dArray& point_array, const double* width, const double* offset, bool relative) {
  spine.cubic_smooth(point_array, relative);
  fill_offsets_and_widths(width, offset);
}

void Path::horizontal(double coord_x, const double* width, const double* offset, bool relative) {
  spine.horizontal(coord_x, relative);
  fill_offsets_and_widths(width, offset);
}
void Path::horizontal(const std::vector<double>& coord_x, const double* width, const double* offset, bool relative) {
  spine.horizontal(coord_x, relative);
  fill_offsets_and_widths(width, offset);
}

void Path::interpolation(const Vec2dArray& point_array, double* angles, bool* angle_constraints, Vec2d* tension, double initial_curl,
                         double final_curl, bool cycle, const double* width, const double* offset, bool relative) {
  spine.interpolation(point_array, angles, angle_constraints, tension, initial_curl, final_curl, cycle, relative);
  fill_offsets_and_widths(width, offset);
}
void Path::mirror(const Vec2d p0, const Vec2d p1) {
  Vec2d v = p1 - p0;
  double tmp = v.length_sq();
  if (tmp == 0)
    return;
  Vec2d r = v * (2 / tmp);
  Vec2d p2 = p0 * 2;
  Vec2d* p = spine.points.data();
  for (uint64_t num = spine.points.size(); num > 0; num--, p++)
    *p = v * (*p - p0).inner(r) - *p + p2;

  for (auto& el : elements) {
    for (auto& wo : el.half_width_and_offset) {
      wo.y = -wo.y;
    }
  }

  // LineElement* el = elements.data();
  // for (uint64_t ne = 0; ne < elements.size(); ne++, el++) {
  //   Vec2d* wo = el->half_width_and_offset.data();
  //   for (uint64_t num = spine.points.size(); num > 0; num--, wo++)
  //     wo->v = -wo->v;
  // }
}
// void Path::offsets(Vec2dArray& result) const {}

void Path::parametric(ParametricVec2 curve_function, void* data, const double* width, const double* offset, bool relative) {
  spine.parametric(curve_function, data, relative);
  fill_offsets_and_widths(width, offset);
}

void Path::quadratic(const Vec2dArray& point_array, const double* width, const double* offset, bool relative) {
  spine.quadratic(point_array, relative);
  fill_offsets_and_widths(width, offset);
}
void Path::quadratic_smooth(Vec2d end_point, const double* width, const double* offset, bool relative) {
  spine.quadratic_smooth(end_point, relative);
  fill_offsets_and_widths(width, offset);
}
void Path::quadratic_smooth(const Vec2dArray& point_array, const double* width, const double* offset, bool relative) {
  spine.quadratic_smooth(point_array, relative);
  fill_offsets_and_widths(width, offset);
}

void Path::rotate(double angle, const Vec2d center) {
  double ca = cos(angle);
  double sa = sin(angle);
  Vec2d* p = spine.points.data();
  for (uint64_t num = spine.points.size(); num > 0; num--, p++) {
    Vec2d q = *p - center;
    p->x = q.x * ca - q.y * sa + center.x;
    p->y = q.x * sa + q.y * ca + center.y;
  }
}

void Path::scale(double scale, const Vec2d center) {
  Vec2d* p = spine.points.data();
  for (uint64_t num = spine.points.size(); num > 0; num--, p++)
    *p = (*p - center) * scale + center;
  Vec2d wo_scale = {1, fabs(scale)};
  if (scale_width)
    wo_scale.u = wo_scale.v;
  LineElement* el = elements.data();
  for (uint64_t ne = 0; ne < elements.size(); ne++, el++) {
    el->end_extensions *= scale;
    Vec2d* wo = el->half_width_and_offset.data();
    for (uint64_t num = spine.points.size(); num > 0; num--)
      *wo++ *= wo_scale;
  }
}

void Path::segment(Vec2d end_point, const double* width, const double* offset, bool relative) {
  spine.segment(end_point, relative);
  // fill_offsets_and_widths(width, offset);
}
void Path::segment(const Vec2dArray& point_array, const double* width, const double* offset, bool relative) {
  spine.segment(point_array, relative);
  // fill_offsets_and_widths(width, offset);
}

gdstk::FlexPath* getFlexPath(Path* path) {
  gdstk::FlexPath* res = (gdstk::FlexPath*)gdstk::allocate_clear(sizeof(gdstk::FlexPath));

  for (auto p : path->spine.points) {
    gdstk::Vec2 pt = {p.x, p.y};
    res->spine.point_array.append(pt);
  }

  res->spine.tolerance = 0.001;
  res->spine.last_ctrl = {0, 0};

  res->num_elements = path->elements.size();
  res->elements = (gdstk::FlexPathElement*)gdstk::allocate_clear(sizeof(gdstk::FlexPathElement) * res->num_elements);

  for (int i = 0; i < path->elements.size(); i++) {
    res->elements[i].half_width_and_offset.ensure_slots(path->elements[i].half_width_and_offset.size());

    for (auto half_and_widht : path->elements[i].half_width_and_offset) {
      gdstk::Vec2 ht = {half_and_widht.x, half_and_widht.y};
      res->elements[i].half_width_and_offset.append(ht);
    }
  }

  res->simple_path = path->simple_path;
  res->scale_width = path->scale_width;

  return res;
}

Polygon::s_ptr Path::to_polygon() {
  PolygonPtrArray polys = {};
  to_polygons(false, qlib::Tag{0, 0}, polys);
  for (auto& p : polys) {
    p->layer = layer;
    p->datatype = datatype;
    return p;
  }
  return nullptr;
}

ErrorCode Path::to_polygons(bool filter, Tag tag, PolygonPtrArray& result) {
  if (width == 0) {
    return ErrorCode::NoError;
  }

  updateElements();

  remove_overlapping_points();
  if (spine.points.size() < 2)
    return ErrorCode::NoError;

  const Array<Vec2d> spine_points = spine.points;
  uint64_t curve_size_guess = spine_points.size() * 2 + 4;

  LineElement* el = elements.data();

  if (elements.size() == 0) {
    qlib::PolygonPtrArray res;
    auto& poly = res.emplace_back();
    poly = qlib::Polygon::create();

    poly->points = spine.points;
    poly->layer = layer;
    poly->datatype = datatype;
    poly->repetition = repetition;
    poly->properties = properties;

    result.push_back(poly);
    return ErrorCode::NoError;
  }

  for (uint64_t ne = 0; ne < elements.size(); ne++, el++) {
    if (filter && (layer != tag.layer || datatype != tag.type))
      continue;

    const double half_widths = el->half_width_and_offset.front().x;
    const double offsets = el->half_width_and_offset.front().y;
    const JoinType join_type = el->join_type;
    const BendType bend_type = el->bend_type;
    const double bend_radius = el->bend_radius;

    Curve right_curve = {};
    Curve left_curve = {};
    right_curve.tolerance = spine.tolerance;
    left_curve.tolerance = spine.tolerance;
    // right_curve.ensure_slots(curve_size_guess);
    // left_curve.ensure_slots(curve_size_guess / 2);

    // Normal to spine segment
    Vec2d spine_normal = (spine_points[1] - spine_points[0]).ortho();
    spine_normal.normalize();
    // First points
    Vec2d p0 = spine_points[0] + spine_normal * offsets;
    Vec2d p1 = spine_points[1] + spine_normal * offsets;
    // Tangent unit vector and segment length
    Vec2d t0 = p1 - p0;
    t0.normalize();
    // Normal to segment
    Vec2d n0 = t0.ortho();

    {  // Initial cap
      const Vec2d cap_l = p0 + n0 * half_widths;
      const Vec2d cap_r = p0 - n0 * half_widths;
      if (el->end_type == EndType::Flush) {
        right_curve.append(cap_l);
        if (half_widths != 0)
          right_curve.append(cap_r);
      } else if (el->end_type == EndType::HalfWidth || el->end_type == EndType::Extended) {
        const double extension = el->end_type == EndType::Extended ? el->end_extensions.u : half_widths;
        if (extension > 0)
          right_curve.append(cap_l);
        right_curve.append(cap_l - extension * t0);
        if (half_widths != 0)
          right_curve.append(cap_r - extension * t0);
        if (extension > 0)
          right_curve.append(cap_r);
      } else if (el->end_type == EndType::Round) {
        right_curve.append(cap_l);
        double initial_angle = n0.angle();
        right_curve.arc(half_widths, half_widths, initial_angle, initial_angle + M_PI, 0);
      } else if (el->end_type == EndType::Smooth) {
        right_curve.append(cap_l);
        const Vec2d p1_l = p1 + n0 * half_widths;
        const Vec2d p1_r = p1 - n0 * half_widths;
        Array<Vec2d> point_array = {};
        point_array.push_back(*(Vec2d*)&cap_r);

        bool angle_constraints[2] = {true, true};
        double angles[2] = {(cap_l - p1_l).angle(), (p1_r - cap_r).angle()};
        Vec2d tension[2] = {Vec2d{1, 1}, Vec2d{1, 1}};
        right_curve.interpolation(point_array, angles, angle_constraints, tension, 1, 1, false, false);
      } else if (el->end_type == EndType::Function) {
        Vec2d dir_l = cap_l - (p1 + n0 * half_widths);
        dir_l.normalize();
        Vec2d dir_r = (p1 - n0 * half_widths) - cap_r;
        dir_r.normalize();
        Array<Vec2d> point_array = (*el->end_function)(cap_l, dir_l, cap_r, dir_r, el->end_function_data);
        right_curve.segment(point_array, false);
        point_array.clear();
      }
    }

    if (spine_points.size() > 2) {
      spine_normal = (spine_points[2] - spine_points[1]).ortho();
      spine_normal.normalize();
      Vec2d p2 = spine_points[1] + spine_normal * offsets;
      Vec2d p3 = spine_points[2] + spine_normal * offsets;
      Vec2d t1 = p3 - p2;
      t1.normalize();
      Vec2d n1 = t1.ortho();
      double u0, u1;
      segments_intersection(p1, t0, p2, t1, u0, u1);
      Vec2d p_next = 0.5 * (p1 + u0 * t0 + p2 + u1 * t1);
      Vec2d p = p0;
      double len_next = (p_next - p).length();

      // Right side: -n
      Vec2d r2 = p - n0 * half_widths;
      Vec2d r3 = p_next - n0 * half_widths;
      Vec2d tr1 = r3 - r2;
      tr1.normalize();

      // Left side: +n
      Vec2d l2 = p + n0 * half_widths;
      Vec2d l3 = p_next + n0 * half_widths;
      Vec2d tl1 = l3 - l2;
      tl1.normalize();

      for (uint64_t i = 1; i < spine_points.size() - 1; i++) {
        Vec2d t2, n2;
        Vec2d r1 = r3;
        Vec2d tr0 = tr1;
        Vec2d l1 = l3;
        Vec2d tl0 = tl1;
        p0 = p2;
        p1 = p3;
        p = p_next;

        if (i + 2 == spine_points.size()) {
          // Last point: no need to find an intersection
          p_next = p1;
          t2 = Vec2d{0, 0};
          n2 = Vec2d{0, 0};
        } else {
          spine_normal = (spine_points[i + 2] - spine_points[i + 1]).ortho();
          spine_normal.normalize();
          p2 = spine_points[i + 1] + spine_normal * offsets;
          p3 = spine_points[i + 2] + spine_normal * offsets;
          t2 = p3 - p2;
          t2.normalize();
          n2 = t2.ortho();
          segments_intersection(p1, t1, p2, t2, u0, u1);
          p_next = 0.5 * (p1 + u0 * t1 + p2 + u1 * t2);
        }

        r2 = p - n1 * half_widths;
        r3 = p_next - n1 * half_widths;
        tr1 = r3 - r2;
        tr1.normalize();

        l2 = p + n1 * half_widths;
        l3 = p_next + n1 * half_widths;
        tl1 = l3 - l2;
        tl1.normalize();

        // Check whether there is enough room for the bend
        double bend_dir = 0;
        double len_factor = 0;
        double center_radius = 0;
        if (bend_type != BendType::None) {
          bend_dir = t0.cross(t1) < 0 ? -1 : 1;
          const Vec2d sum_t = t0 + t1;
          const double len_prev = len_next;
          len_next = (p_next - p).length();
          len_factor = (fabs(sum_t.x) > fabs(sum_t.y)) ? bend_dir * (n0.x - n1.x) / sum_t.x : bend_dir * (n0.y - n1.y) / sum_t.y;
          center_radius = bend_radius - bend_dir * offsets;
          const double len_required = len_factor * center_radius;
          if (len_required > len_prev || len_required > len_next || center_radius <= half_widths) {
            // Not enough room for the bend
            bend_dir = 0;
          } else {
            len_next -= len_required;
          }
        }

        if (bend_dir < 0) {
          const double initial_angle = n0.angle();
          double final_angle = n1.angle();
          if (final_angle > initial_angle)
            final_angle -= 2 * M_PI;
          const Vec2d center = p - 0.5 * center_radius * (n0 + n1 + len_factor * (t0 - t1));

          // Right: inner side of the bend
          double radius = center_radius - half_widths;
          if (bend_type == BendType::Circular) {
            const Vec2d arc_start = center + n0 * radius;
            right_curve.append(arc_start);
            right_curve.arc(radius, radius, initial_angle, final_angle, 0);
          } else if (bend_type == BendType::Function) {
            Array<Vec2d> point_array = (*el->bend_function)(radius, initial_angle, final_angle, center, el->bend_function_data);
            right_curve.segment(point_array, false);
            point_array.clear();
          }

          // Left: outer side of the bend
          radius = center_radius + half_widths;
          if (bend_type == BendType::Circular) {
            const Vec2d arc_start = center + n0 * radius;
            left_curve.append(arc_start);
            left_curve.arc(radius, radius, initial_angle, final_angle, 0);
          } else if (bend_type == BendType::Function) {
            Array<Vec2d> point_array = (*el->bend_function)(radius, initial_angle, final_angle, center, el->bend_function_data);
            left_curve.segment(point_array, false);
            point_array.clear();
          }
        } else if (bend_dir > 0) {
          const double initial_angle = (-n0).angle();
          double final_angle = (-n1).angle();
          if (final_angle < initial_angle)
            final_angle += 2 * M_PI;
          Vec2d center = p + 0.5 * center_radius * (n0 + n1 + len_factor * (t1 - t0));

          // Right: outer side of the bend
          double radius = center_radius + half_widths;
          if (bend_type == BendType::Circular) {
            const Vec2d arc_start = center - n0 * radius;
            right_curve.append(arc_start);
            right_curve.arc(radius, radius, initial_angle, final_angle, 0);
          } else if (bend_type == BendType::Function) {
            Array<Vec2d> point_array = (*el->bend_function)(radius, initial_angle, final_angle, center, el->bend_function_data);
            right_curve.segment(point_array, false);
            point_array.clear();
          }

          // Left: inner side of the bend
          radius = center_radius - half_widths;
          if (bend_type == BendType::Circular) {
            const Vec2d arc_start = center - n0 * radius;
            left_curve.append(arc_start);
            left_curve.arc(radius, radius, initial_angle, final_angle, 0);
          } else if (bend_type == BendType::Function) {
            Array<Vec2d> point_array = (*el->bend_function)(radius, initial_angle, final_angle, center, el->bend_function_data);
            left_curve.segment(point_array, false);
            point_array.clear();
          }
        } else {
          if (tr0.cross(tr1) < 0) {
            // Right: inner side of the bend
            segments_intersection(r1, tr0, r2, tr1, u0, u1);
            const Vec2d ri = 0.5 * (r1 + u0 * tr0 + r2 + u1 * tr1);
            right_curve.append(ri);
          } else {
            // Right: outer side of the bend
            if (join_type == JoinType::Bevel) {
              right_curve.append(r1);
              right_curve.append(r2);
            } else if (join_type == JoinType::Miter) {
              segments_intersection(r1, tr0, r2, tr1, u0, u1);
              const Vec2d ri = 0.5 * (r1 + u0 * tr0 + r2 + u1 * tr1);
              right_curve.append(ri);
            } else if (join_type == JoinType::Natural) {
              segments_intersection(r1, tr0, r2, tr1, u0, u1);
              const double half_width = half_widths;
              u1 = -u1;
              if (u0 <= half_width && u1 <= half_width) {
                const Vec2d ri = 0.5 * (r1 + u0 * tr0 + r2 - u1 * tr1);
                right_curve.append(ri);
              } else {
                const Vec2d ri0 = r1 + (u0 > half_width ? half_width : u0) * tr0;
                right_curve.append(ri0);
                const Vec2d ri1 = r2 - (u1 > half_width ? half_width : u1) * tr1;
                right_curve.append(ri1);
              }
            } else if (join_type == JoinType::Round) {
              right_curve.append(r1);
              const double initial_angle = (-n0).angle();
              double final_angle = (-n1).angle();
              if (final_angle < initial_angle)
                final_angle += 2 * M_PI;
              right_curve.arc(half_widths, half_widths, initial_angle, final_angle, 0);
            } else if (join_type == JoinType::Smooth) {
              right_curve.append(r1);
              Array<Vec2d> point_array = {};
              point_array.push_back(*(Vec2d*)&r2);
              bool angle_constraints[2] = {true, true};
              double angles[2] = {tr0.angle(), tr1.angle()};
              Vec2d tension[2] = {Vec2d{1, 1}, Vec2d{1, 1}};
              right_curve.interpolation(point_array, angles, angle_constraints, tension, 1, 1, false, false);
            } else if (join_type == JoinType::Function) {
              Array<Vec2d> point_array = (*el->join_function)(r1, tr0, r2, tr1, p, half_widths * 2, el->join_function_data);
              right_curve.segment(point_array, false);
              point_array.clear();
            }
          }

          if (tl0.cross(tl1) > 0) {
            // Left: inner side of the bend
            segments_intersection(l1, tl0, l2, tl1, u0, u1);
            const Vec2d li = 0.5 * (l1 + u0 * tl0 + l2 + u1 * tl1);
            left_curve.append(li);
          } else {
            // Left: outer side of the bend
            if (join_type == JoinType::Bevel) {
              left_curve.append(l1);
              left_curve.append(l2);
            } else if (join_type == JoinType::Miter) {
              segments_intersection(l1, tl0, l2, tl1, u0, u1);
              const Vec2d li = 0.5 * (l1 + u0 * tl0 + l2 + u1 * tl1);
              left_curve.append(li);
            } else if (join_type == JoinType::Natural) {
              segments_intersection(l1, tl0, l2, tl1, u0, u1);
              const double half_width = half_widths;
              u1 = -u1;
              if (u0 <= half_width && u1 <= half_width) {
                const Vec2d li = 0.5 * (l1 + u0 * tl0 + l2 - u1 * tl1);
                left_curve.append(li);
              } else {
                const Vec2d li0 = l1 + (u0 > half_width ? half_width : u0) * tl0;
                left_curve.append(li0);
                const Vec2d li1 = l2 - (u1 > half_width ? half_width : u1) * tl1;
                left_curve.append(li1);
              }
            } else if (join_type == JoinType::Round) {
              left_curve.append(l1);
              const double initial_angle = n0.angle();
              double final_angle = n1.angle();
              if (final_angle > initial_angle)
                final_angle -= 2 * M_PI;
              left_curve.arc(half_widths, half_widths, initial_angle, final_angle, 0);
            } else if (join_type == JoinType::Smooth) {
              left_curve.append(l1);
              Array<Vec2d> point_array = {};
              point_array.push_back(*(Vec2d*)&l2);
              bool angle_constraints[2] = {true, true};
              double angles[2] = {tl0.angle(), tl1.angle()};
              Vec2d tension[2] = {Vec2d{1, 1}, Vec2d{1, 1}};
              left_curve.interpolation(point_array, angles, angle_constraints, tension, 1, 1, false, false);
            } else if (join_type == JoinType::Function) {
              Array<Vec2d> point_array = (*el->join_function)(l1, tl0, l2, tl1, p, half_widths * 2, el->join_function_data);
              left_curve.segment(point_array, false);
              point_array.clear();
            }
          }
        }

        t0 = t1;
        n0 = n1;
        t1 = t2;
        n1 = n2;
      }
    }

    {  // End cap
      const uint64_t last = spine_points.size() - 1;
      const Vec2d cap_l = p1 + n0 * half_widths;
      const Vec2d cap_r = p1 - n0 * half_widths;
      if (el->end_type == EndType::Flush) {
        left_curve.append(cap_l);
        if (half_widths != 0)
          left_curve.append(cap_r);
      } else if (el->end_type == EndType::HalfWidth || el->end_type == EndType::Extended) {
        const double extension = el->end_type == EndType::Extended ? el->end_extensions.v : half_widths;
        if (extension > 0)
          left_curve.append(cap_l);
        left_curve.append(cap_l + extension * t0);
        if (half_widths != 0)
          left_curve.append(cap_r + extension * t0);
        if (extension > 0)
          left_curve.append(cap_r);
      } else if (el->end_type == EndType::Round) {
        left_curve.append(cap_l);
        double initial_angle = n0.angle();
        left_curve.arc(half_widths, half_widths, initial_angle, initial_angle - M_PI, 0);
      } else if (el->end_type == EndType::Smooth) {
        left_curve.append(cap_l);
        const Vec2d p0_l = p0 + n0 * half_widths;
        const Vec2d p0_r = p0 - n0 * half_widths;
        Array<Vec2d> point_array = {};
        point_array.push_back(*(Vec2d*)&cap_r);
        bool angle_constraints[2] = {true, true};
        double angles[2] = {(cap_l - p0_l).angle(), (p0_r - cap_r).angle()};
        Vec2d tension[2] = {Vec2d{1, 1}, Vec2d{1, 1}};
        left_curve.interpolation(point_array, angles, angle_constraints, tension, 1, 1, false, false);
      } else if (el->end_type == EndType::Function) {
        Vec2d dir_r = cap_r - (p0 - n0 * half_widths);
        dir_r.normalize();
        Vec2d dir_l = (p0 + n0 * half_widths) - cap_l;
        dir_l.normalize();
        Array<Vec2d> point_array = (*el->end_function)(cap_r, dir_r, cap_l, dir_l, el->end_function_data);
        const uint64_t count = point_array.size();
        for (uint64_t j = 0; j < count / 2; j++) {
          Vec2d tmp = point_array[count - 1 - j];
          point_array[count - 1 - j] = point_array[j];
          point_array[j] = tmp;
        }
        left_curve.segment(point_array, false);
        point_array.clear();
      }
    }

    int old_size = right_curve.ensure_slots(left_curve.points.size());
    Vec2d* dst = right_curve.points.data() + old_size;
    Vec2d* src = left_curve.points.data() + left_curve.points.size() - 1;
    for (uint64_t i = left_curve.points.size(); i > 0; i--)
      *dst++ = *src--;

    left_curve = Curve();

    curve_size_guess = right_curve.points.size() * 6 / 5;

    // 过滤计算过程中产生的无穷点位信息
    qlib::Vec2dArray new_points;
    if (!right_curve.points.empty()) {
      for (auto p : right_curve.points) {
        if (!std::isinf(p.x) && !std::isinf(p.y) && !std::isnan(p.x) && !std::isnan(p.y)) {
          new_points.push_back(p);
        }
      }
    }
    Polygon::s_ptr result_polygon = Polygon::create();
    result_polygon->layer = layer;
    result_polygon->datatype = datatype;
    result_polygon->points = new_points;
    result_polygon->repetition = repetition;
    result_polygon->properties = properties;

    result.push_back(result_polygon);
  }
  return ErrorCode::NoError;
}

void Path::translate(double dx, double dy) {
  Vec2d* p = spine.points.data();
  Vec2d d = Vec2d(dx, dy);
  for (uint64_t num = spine.points.size(); num > 0; num--)
    *p++ += d;
}

void Path::turn(double radius, double angle, const double* width, const double* offset) {
  spine.turn(radius, angle);
  fill_offsets_and_widths(width, offset);
}
void Path::vertical(double coord_y, const double* width, const double* offset, bool relative) {
  spine.vertical(coord_y, relative);
  fill_offsets_and_widths(width, offset);
}
void Path::vertical(const std::vector<double>& coord_y, const double* width, const double* offset, bool relative) {
  spine.vertical(coord_y, relative);
  fill_offsets_and_widths(width, offset);
}

void Path::widths(std::vector<double>& res) {}

void Path::updateType() {
  if (elements.size() == 1 && elements[0].half_width_and_offset[0].e[0] == 0) {
    type = PathType::Spine;
  }
  if (elements.size() == 1 && elements[0].half_width_and_offset[0].e[0] > 0) {
    type = PathType::One;
  }
  if (elements.size() == 2 && elements[0].half_width_and_offset[0].e[0] > 0) {
    type = PathType::WaveGuide;
  }

  auto p = Property::create();
  auto p_v = PropertyValue::create();
  p->values.push_back(p_v);
  p->name = "type";
  p_v->type = PropertyType::UnsignedInteger;
  p_v->unsigned_integer = type;

  if (properties.size() == 0) {
    properties.push_back(p);
  } else {
    properties[0] = p;
  }
}

void Path::updateElements() {
  // 拉伸线只能拉伸 无宽度无中心线宽的线
  // 测试要求修改共面波导宽度和中心线宽为0后可以拉伸
  drawType();

  elements.resize(2);
  {
    elements[0].bend_radius = radius;
    if (elements[0].bend_radius > 0) {
      elements[0].bend_type = BendType::Circular;
      elements[0].join_type = JoinType::Round;
    } else {
      elements[0].bend_type = BendType::None;
      elements[0].join_type = JoinType::Natural;
    }

    elements[0].half_width_and_offset.resize(1);
    elements[0].half_width_and_offset[0] = Vec2d(width / 2.0, width / 2.0 + center_width / 2.0);
  }

  {
    elements[1].bend_radius = radius;
    if (elements[1].bend_radius > 0) {
      elements[1].bend_type = BendType::Circular;
      elements[1].join_type = JoinType::Round;
    } else {
      elements[1].bend_type = BendType::None;
      elements[1].join_type = JoinType::Natural;
    }

    elements[1].half_width_and_offset.resize(1);
    elements[1].half_width_and_offset[0] = Vec2d(width / 2.0, -(width / 2.0 + center_width / 2.0));
  }
}

void Path::update() {
  updateElements();
  // updateType();
}

// 计算包围盒(min,max)
Vec2dArray Path::bounding_box() {
  Vec2dArray res(2);

  Vec2d& min = res[0];
  Vec2d& max = res[1];
  min.x = min.y = DBL_MAX;
  max.x = max.y = -DBL_MAX;

  auto points = spine.points;
  if (width > 0) {
    PolygonPtrArray polys = {};
    to_polygons(false, Tag{0, 0}, polys);
    for (auto& p : polys) {
      points.insert(points.end(), p->points.begin(), p->points.end());
    }
  }

  // 计算未应用重复下的包围盒
  // auto& points = spine.points;
  const Vec2d* p = points.data();
  for (uint64_t num = points.size(); num > 0; num--, p++) {
    if (p->x < min.x)
      min.x = p->x;
    if (p->x > max.x)
      max.x = p->x;
    if (p->y < min.y)
      min.y = p->y;
    if (p->y > max.y)
      max.y = p->y;
  }

  // 考虑重复
  if (repetition.type != RepetitionType::None) {
    // 获取重复下的极值
    Vec2dArray offsets = {};
    repetition.get_extrema(offsets);

    Vec2d* off = offsets.data();
    Vec2d min0 = min;
    Vec2d max0 = max;

    // 将极值应用到每个包围盒上
    for (uint64_t i = offsets.size(); i > 0; i--, off++) {
      if (min0.x + off->x < min.x)
        min.x = min0.x + off->x;
      if (max0.x + off->x > max.x)
        max.x = max0.x + off->x;
      if (min0.y + off->y < min.y)
        min.y = min0.y + off->y;
      if (max0.y + off->y > max.y)
        max.y = max0.y + off->y;
    }
  }

  return res;
}

void Path::transform(double magnification, bool x_reflection, double rotation, const Vec2d origin) {
  double ca = cos(rotation);
  double sa = sin(rotation);
  for (auto& p : spine.points) {
    Vec2d q = p * magnification;
    if (x_reflection)
      q.y = -q.y;
    p.x = q.x * ca - q.y * sa + origin.x;
    p.y = q.x * sa + q.y * ca + origin.y;
  }
}

Vec2dArray Path::render_aabb() {
  return bounding_box();
}

Vec2dArray Path::render_poly() {
  return spine.points;
}

Vec2dArray Path::render_absorp() {
  return spine.points;
}

void Path::updateWidthAndRadius() {
  elements.resize(1);
  elements[0].bend_radius = radius;

  if (elements[0].bend_radius > 0) {
    elements[0].bend_type = BendType::Circular;
    elements[0].join_type = JoinType::Round;
  } else {
    elements[0].bend_type = BendType::None;
    elements[0].join_type = JoinType::Natural;
  }

  // elements[0].end_type = EndType::HalfWidth;
  // elements[0].join_type = JoinType::Bevel;
  // elements[0].bend_type = qlib::BendType::Circular;

  elements[0].half_width_and_offset.resize(1);

  for (size_t i = 0; i < 1; i++) {
    elements[0].half_width_and_offset[i] = Vec2d(width / 2.0, 0);
  }
}

void Path::widthToElement() {
  elements.resize(1);
  elements[0].bend_radius = radius;

  if (elements[0].bend_radius > 0) {
    elements[0].bend_type = BendType::Circular;
    elements[0].join_type = JoinType::Round;
  } else {
    elements[0].bend_type = BendType::None;
    elements[0].join_type = JoinType::Natural;
  }

  // elements[0].end_type = EndType::HalfWidth;
  // elements[0].join_type = JoinType::Bevel;
  // elements[0].bend_type = qlib::BendType::Circular;

  elements[0].half_width_and_offset.resize(1);

  for (size_t i = 0; i < 1; i++) {
    elements[0].half_width_and_offset[i] = Vec2d(width / 2.0, 0);
  }
}

double Path::getLineWidth() {
  return width;
}
void Path::setLineWidth(double _width) {
  width = _width;
  update();
}

double Path::getLineRadius() {
  return radius;
}
void Path::setLineRadius(double _radius) {
  radius = _radius;
  update();
}

uint32_t Path::getLayer() {
  return layer;
}
void Path::setLayer(uint32_t _layer) {
  layer = _layer;
}

PathParts Path::get_parts() {
  auto v_points = spine.points;
  auto v_edges = get_edges(v_points);
  auto v_edges_mid = get_edges_mid(v_edges);

  auto box = bounding_box();
  auto box_points = box_to_points(box);
  auto box_edges = get_edges(box_points);
  auto box_edges_mid = get_edges_mid(box_edges);

  auto center = box_center(box);
  auto res = std::make_tuple(center, v_points, v_edges, v_edges_mid, box_points, box_edges, box_edges_mid);
  return res;
}

Vec2dArray Path::query_parts(GeoPartType type, uint32_t id) {
  auto parts = get_parts();
  if (type == GeoPartType::Center)
    return Vec2dArray{std::get<0>(parts)};
  if (type == GeoPartType::Point)
    return Vec2dArray{std::get<1>(parts)[id]};
  if (type == GeoPartType::Edge)
    return std::get<2>(parts)[id];
  if (type == GeoPartType::EdgeMid)
    return Vec2dArray{std::get<3>(parts)[id]};
  if (type == GeoPartType::BoxPoint)
    return Vec2dArray{std::get<4>(parts)[id]};
  if (type == GeoPartType::BoxEdge)
    return std::get<5>(parts)[id];
  if (type == GeoPartType::BoxEdgeMid)
    return Vec2dArray{std::get<6>(parts)[id]};
  return Vec2dArray{};
}

void Path::coplanarWaveguide(double line_width, double center_line_width, double bend_radius, uint32_t layer_) {
  width = line_width;
  center_width = center_line_width;
  radius = bend_radius;
  layer = layer_;
  type = PathType::WaveGuide;
  update();
}

double Path::getCenterWidth() {
  return center_width;
}

void Path::setCenterWidth(double width) {
  center_width = width;
  update();
}

void Path::drawType() {
  if (width == 0.0) {
    type = PathType::Spine;
    return;
  }

  if (center_width != 0.0) {
    type = PathType::WaveGuide;
  } else if (width != 0.0) {
    type = PathType::One;
  } else {
    type = PathType::Spine;
  }
}

void Path::setNetName(std::string net_name) {
  m_net_name = net_name;
  set_gds_property(this->properties, 1000, net_name);
}

std::string Path::getNetName() {
  auto values = get_gds_property(this->properties, 1000);
  if (values) {
    for (auto v : *values) {
      if (v->string != "") {
        return v->string;
      }
    }
  }
  return std::string();
}
}  // namespace qlib